Numerical simulation evidence of spectrum rearrangement in impure graphene 
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By means of numerical simulation we confirm that in graphene with point defects a quasigap 
opens in the vicinity of the resonance state with increasing impurity concentration. We prove that 
states inside this quasigap cannot longer be described by a wavevector and are strongly localized. 
We visualize states corresponding to the density of states maxima within the quasigap and show 
that they are yielded by impurity pair clusters. 
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I. INTRODUCTION 



Not so long ago, graphene has been cleaved out for 
the first time by the so-called scotch-tape technique.- 
Since that, this novel, truly two-dimensional material 
continues to put new challenges to scientific community. 
Graphene is already known to manifest some remarkable 
properties. The most unusual of them, and, correspond- 
ingly, the most frequently emphasized on, is the Dirac 
dispersion of Fermi elementary excitations. This unusual 
spectrum makes graphene a rather promising material 
for a variety of applications in computer electronics and 
chemical sensors. While graphene is known to possess 
outstanding structural stability and crystalline quality, 
existing methods of its isolation necessarily imply the 
presence of a certain amount of defects. Moreover, im- 
purities can be embedded into graphene intentionally in 
order to tune up its physical characteristics in accordance 
with a specific application. Even though some applica- 
tions are destined for a distant future, the need in delib- 
erate and proper functionalization of graphene provides 
adequate grounds for an extensive study of different types 
of defects in this material. Despite a noticeable quantity 
of papers devoted to the study of impure graphene, a 
comprehensive understanding of impurity effects on its 
electron spectrum is still lacking. 

While different types of disorder are inherent in 
graphene, below we are particularly interested among 
them in the substitutional point defects. This commonly 
used model applies not only to chemically substituted 
carbon atoms or an absence of them (vacancies), but 
also, to a known extent, to adsorbed atoms, molecules, 
or radicals on the graphene sheet. 2 Regarding the substi- 
tutional impurities in graphene, comprehensive attention 
has been paid as to the single impurity problem, in which 
the impurity state wave function has been studied for a 
single defect and a pair of them, 3 as to the evolution of 
the density of states (DOS) in graphene with increasing 
the impurity concentration4^^^^ However, such a phe- 
nomenon as spectrum rearrangement is frequently over- 
looked. The main concept of the spectrum rearrangement 



is based on the existence of some critical impurity con- 
centration, at which the spectrum of a disordered system 
undergoes a cardinal qualitative change. As a rule, the 
spectrum rearrangement should be related to the appear- 
ance of a local level or a resonance state. Characteristics 
of this impurity state, namely its energy and damping, 
which are determined within the single impurity prob- 
lem, are shaping the scenario and type of the subsequent 
spectrum rearrangement. Albeit a single point defect is 
expected to perturb only the lattice cite occupied by an 
impurity or, at most, the adjacent lattice cites, and thus 
should be classified as a short-range defect, the effective 
radius of the correspondent impurity state can far exceed 
the lattice constant, when its energy is close to the van 
Hove singularity in the spectrum. Spacial overlap of indi- 
vidual impurity states, which occurs with increasing the 
impurity concentration, is marking the onset of the spec- 
trum rearrangement. This simple consideration gives a 
possibility to roughly estimate the critical concentration 
of the spectrum rearrangement. Since the effective radius 
of the single impurity state in certain cases can be large 
compared to the lattice constant, the respective critical 
concentration should be fairly low. As a result, in such 
situations only a trace amount of impurities can provoke 
the spectrum rearrangement. 

In graphene, the spectrum rearrangement driven by 
an increase in the concentration of defects described by 
the Lifshitz model has been analytically examined in 
Refs. Isljlfj. It has been demonstrated that the passage 
of the spectrum rearrangement is of the anomalous type 
due to a weak resonance state. That means that for 
a sufficiently large compared to the bandwidth impurity 
potential a quasigap is gradually opening around the res- 
onance energy and at the critical concentration spans up 
to the Dirac point in the spectrum. With further in- 
crease in the impurity concentration (i.e. after the spec- 
trum rearrangement), the quasigap is rapidly enlarging. 
In this regime the width of the quasigap is approximately 
proportional to the square root of the impurity concen- 
tration. In Refs. |Q|J 10 the analysis of the spectrum rear- 
rangement in graphene has been conducted by means of 
the coherent potential approximation (CPA) applicabil- 
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ity criterion, which is instrumental in determining spec- 
trum domains with different degree of localization^ 

The aim of the current work is to carry out the numer- 
ical simulation of graphene DOS with the special em- 
phasis on the spectrum rearrangement phenomenon (for 
its physical description see also the review Ref . [l2l) . By 
considering for each chosen perturbation strength those 
impurity concentrations that are close to the expected 
critical concentration of the spectrum rearrangement, we 
show that implementing the criterion of the CPA appli- 
cability it is possible to judge upon the spectrum rear- 
rangement process in graphene. As a next step in our 
previous attempts^*^ we are paying a special attention 
to the CPA validity. By comparing the CPA output with 
the numerical results, we verify that the CPA applica- 
bility criterion works rather well. We also notice that 
when impurity concentration is low enough, the aver- 
age T-matrix approximation (ATA) is in good agreement 
with numerically calculated DOS. We discuss the spec- 
trum rearrangement in graphene and the correspondent 
interplay between numerical and analytical results. After 
the spectrum rearrangement takes place, we identify the 
cluster structure of graphene's DOS in the vicinity of the 
impurity resonance energy. The structure observed evi- 
dently cannot be explained by means of available analytic 
approaches and requires further analysis. 

As an effective tool for numerical calculations we em- 
ploy the negative eigenvalue theorem as suggested by 
Dean»i2 This approach, to the best of our knowledge, 
hasn't been used for graphene yet (see, e.g. Refs. 
Its advantages are discussed below. Finally, we describe 
a drift of the Fermi energy from the DOS minimum, a 
kind of a self-doping effect, when Fermi level shifts away 
without actual introduction of additional carriers into the 
disordered system. 

The paper is organized as follows: in Section II A we 
remind the basic mathematical impurity model. In Sec- 
tion II B we introduce the concept of spectrum rearrange- 
ment. In Section II C we briefly set forth the numerical 
approach. In Section III we present and discuss results. 
In the last Section we summarize outcomes of our study. 



II. MODEL DISORDERED SYSTEM AND 
METHODS OF ITS ANALYSIS 

A. Impurity model 

In the tight-binding approximation the simplest (for 
spinless fermions) graphene Hamiltonian has the form, 
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Here, t rs 2.7eV is the hopping integral, n and m denote 
vectors of lattice cells, Greek indices a and (3 correspond 
to the graphene sublattices, and summation runs over all 
nearest neighbors. It has been confirmed experimentally 
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FIG. 1: (color online) A comparison between the exact go(e) 
(solid line) and its analytical approximation © (dash-dotted 
line) . 



that this approximation describes graphene's electronic 
spectrum fairly wellJ^ 

The diagonal element of the Green's function in the 
cite representation, 



(e) = lim 



(5->o+ s + id — e(j) 



(2) 



where the summation runs over all eigenstates and 
e(j) are their corresponding eigenvalues, in the case of 
Hamiltonian |l} can be easily approximated by 



ffo(e) = .9o Q o Q (e) 



W 2 \W 



w 2 



(3) 



in the low-energy limit, i.e. |e| <C W = y 7r\/3t, where 
W is the bandwidth. A detailed derivation of ([3]) can be 
found in Refs. l9lficl A comparison between the exact di- 
agonal element of the Green's function and its low-energy 
approximation is given in Fig. [T] Since we are interested 
only in the relatively close vicinity of the Dirac point in 
the spectrum, the approximation (|3|) looks appropriately 
shaped. 

While the host DOS can be straightforwardly found 
from the imaginary part of the latter expression, impu- 
rities break the translational symmetry and so the DOS 
of disordered graphene (which is the focus of the current 
investigation) cannot be directly obtained. Point defects 
in graphene are usually modeled by adding the following 
perturbation in the Hamiltonian (the so-called Lifshitz 
model): 



(4) 



where Vl is the impurity potential, and p a runs over sites 
occupied by impurities. It is supposed that impurities 
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are distributed among lattice sites absolutely at random 
with concentration c representing the probability that 
an arbitrary site is occupied by an impurity. Thus, for 
a large system with N lattice sites the total amount of 
impurities tends to cN. 



B. Spectrum rearrangement and CPA applicability 
criterion 

When the impurity concentration is small enough, con- 
ventional analytic approache d 15 ' 16 can be applied. To 
be concise, in the first approximation the averaged per- 
turbed Green's function 



leads to the intensification of the impurity state local- 



G(e) =< (e-Ho-U)- 1 > 



(5) 



can be found as some renormalization of the host Green's 
function: 



G(e) = g( e - £ r(e)) 

Two cases are of a particular interest: 
matrix approximation (ATA), 



(6) 

the average T- 



<jata(£) = 



cV L 



I - (1 - c)V L g (s) 
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which accounts for multiple single-site scattering by an 
impurity, and the coherent potential approximation, 



ctcpa(£) = 



cV L 



1 - [Vl — (TcpA{s)]go{e - ucpa{£)) ' 



(8) 



which adds the self-consistency. In the CPA the self- 
energy is taken from the requirement that the single-site 
renormalized T-matrix should be zero on average. In 
both methods scatterings on pairs and larger groups of 
impurities is omitted. Thus, these approximations are 
expected to remain valid, when cluster effects are insignif- 
icant in a disordered system. 

However, when impurity concentration is gradually 
increased, individual impurity states (visualized for 
graphene, e.g., in Ref. 3) begin to overlap with each 
other. Thus, a contribution from scatterings on impu- 
rity clusters to the self-energy is becoming more pro- 
nounced in the vicinity of the impurity state energy. As 
a result, a significant overlap between impurity states 
corresponds to the commencement of substantial mod- 
ifications in the spectrum of a disordered system. In 
other words, it points out the critical concentration of 
the spectrum rearrangement. This simple reasoning pro- 
vides a possibility to estimate the critical concentration 
in graphene with Lifshitz impurities. From the expres- 
sion for the non-diagonal element of the host Green's 
function,^ it can be deduced that an effective decay ra- 



dius of the impurity state is r,. 



\V L \. It should be 



noticed that in commonly encountered cases for the Lif- 
shitz impurity model, an increase in the parameter Vl 



ization and, consequently, to a decrease in 



imp • 



The 



opposite result for graphene is caused by the particle- 
hole symmetry of the Dirac Hamiltonian. Another char- 
acteristic space interval is the average distance between 
impurities, which depends on impurity concentration as 



< r >~ \j\fc. Both radii coincide (< r >? 



,) at 



some impurity concentration c r ~ Thus, the con- 

dition < r >« r imp defines the spectrum rearrangement 
concentration c r . As has been argued in Ref. [^, at this 
critical concentration a quasigap filled with strongly lo- 
calized states should sweep from the resonance energy 
to the Dirac point, stimulating an accumulation of con- 
siderable changes in the DOS. Albeit this condition is 
reasonably intuitive, it is too rough for an accurate fore- 
cast of the critical concentration. Moreover, it does not 
provide any information on the energy intervals, in which 
the spectrum is more exposed to the rearrangement pro- 
cess. Since cluster effects, which make up the core of 
the spectrum rearrangement, are not included into the 
CPA, the CPA applicability criterion can be successfully 
employed for the analysis of the spectrum rearrangement 
process. This very approach has been in fact realized in 
Refs.lMIolfor the impure graphene. Namely, by following 
the conventional technique of the Green's function cluster 
expansio n 12 ' 16 it is possible to represent the self-energy 
as a series in all possible groups of impurity centers. At 
that, the first term of this series corresponds to the con- 
ventional CPA. The small parameter of the series, 



R(e) 



V L - (tcp a (e) 
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— o"cpa(£) 



1 + ocpA(e)ffo(£ ~ a CPA{£)) 
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is indicative of the relative strength of cluster effects at 
a given energy, and can be used to outline qualitatively 
different spectrum domains. 

In those spectral domains, where the small parameter 
of the series R(e) is high, the CPA is not reliable and cor- 
respondent states are showing a tendency towards local- 
ization. For ordinary 3D systems a mobility edge should 
be expected close to the energy, at which R(e) — 
At the same time, it can be demonstrated that the maxi- 
mum magnitude of i?(e) is unity, and it is reached at the 
van Hove singularities of the spectrum. At the low im- 
purity concentration, expression ^ can be approximated 
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FIG. 2: (color online) Density of states for graphene without 
impurities obtained by Dean's numerical method compared 
to the exact one for the infinite system. 



by 



R(e) = 



'CPA 



9o(E) 



dgo(E) 



dE 
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There are two factors in ([TO]) . The first of them, 
a CPA( £ )/ c i increases in absolute magnitude around the 
impurity state energy, which can be determined from the 
Lifshitz equation, 1 = Vl.9o(£))- The second one (in the 
square brackets) is the sum of the Green's function and 
its derivative, which increases in the vicinity of the Dirac 
point (or any other van Hove singularity) . Consequently, 
the energy dependence of the CPA applicability criterion 
should possess different maxima, around which the CPA 
is not valid. Even though the CPA applicability criterion 
has been deduced for the fictitious system with a single 
Dirac cone in the spectrum, it will be apparent below 
that it is an adequate tool for the spectrum rearrange- 
ment analysis in the actual graphene. As regards the 
CPA and the ATA, it is not difficult to show that pres- 
ence of the two different Dirac cones in graphene alters 
their output only in a trivial way. 



C. Numerical method 

Numerical techniques involving diagonalization of the 
random matrix are too resource consuming to simulate 
disordered systems approaching in their dimensions real 
experimental samples. However, information on eigen- 
vectors is superficial for the DOS calculations. So far, 
the Haydock method^ based on an expansion of the di- 
agonal element of the Green's function into an infinite 
fraction, has been extensively used for numerical calcu- 
lation of the graphene DOS*££ Still, this approach is not 
without its shortcomings. It is the local DOS (LDOS) 
that is calculated within this approach. The necessity 
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FIG. 3: (color online) A set of figures corresponding to impu- 
rity perturbation Vl = 4t and different concentrations. Crit- 
ical concentration is c r « 0.015. Stepped curve stands for 
the numerical computation, dashed — the CPA, dash-dotted 
— the ATA (left Y-axis represents their values). Solid black 
curve is R(e) (right Y-axis represents its values) . Triangle on 
the energy axis denotes the Fermi level position. 
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FIG. 4: (color online) A set of figures corresponding to impurity perturbation Vl = St and different concentrations. Critical 
concentration is c r ~ 0.003. Stepped curve stands for the numerical computation, Dashed — the CPA, dash-dotted — the 
ATA (left Y-axis represents their values). Solid black curve is R(e) (right Y-axis represents its values). Triangle on the energy 
axis denotes the Fermi level position. 



to truncate the infinite fraction at some point sometimes 
leads to unphysical oscillations of the LDOS, which are 
difficult to keep under control and to distinguish from 
actual features of the spectrum. Furthermore, the total 
DOS is obtained in the Haydock method by averaging 
the LDOS at several lattice sites, and an inclusion of all 
sites in the model system into the averaging process is 
absolutely impractical. The above leaves a touch of un- 
certainty in the DOS minutiae. In a contrast, we relied on 
the Dean's calculation schemed It allows to obtain the 
total number of eigenvalues of a hermitian matrix that 
are less than a specified value. This provides a possibil- 
ity to explore the DOS with a desired degree of precision 
and to preserve all particularities of the resulting curve. 

This method has proven to be especially effective for 
ID systems. The time required to finish a single Dean's 
algorithm loop is proportional to the number of atoms 
(N) in a ID system, which is fast enough to simulate 
really large ID systems. However, with an increase in 



the system's dimensionality, the computational time re- 
quired for one loop increases. For a 2D system it is pro- 
portional to N 2 M- In our case of graphene, we obtained 
DOSs for the system comprised of 5.3 • 10 6 atoms, which 
corresponds to a system with the linear dimensions about 
0.3/im — about the size of real experimental samples. To 
eliminate the influence of boundary states on the DOS 
we applied periodical boundary condition for the zigzag 
boundary of the model system under consideration^ 
The numerically calculated DOS for the described model 
system is given in Fig. [2] Some jaggedness seen in the 
DOS curve is related to the finite size of the model sys- 
tem. 



III. RESULTS AND DISCUSSION 

The figures Fig. [3J Fig. 0] and Fig. [5] correspond to 
impurity perturbations Vl = it, Vl = St, and Vl = 16i , 
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FIG. 5: (color online) A set of figures corresponding to the impurity perturbation Vl = 16t and different concentrations. 
Critical concentration is c r ~ 0.0007. Stepped curve stands for the numerical computation, dashed — the CPA, dash-dotted 
— the ATA (left Y-axis represents their values) . Solid black curve is R(e) (right Y-axis represents its values) . Triangle on the 
energy axis denotes the Fermi level position. 



respectively. At the negative impurity potentials Vl the 
whole picture is simply mirrored against the zero energy. 
For each perturbation magnitude we consider qualita- 
tively different regions of impurity concentration: c < c r 
(before the spectrum rearrangement), c ss c r (in the 
course of the spectrum rearrangement) and c > c r (af- 
ter the spectrum rearrangement). We plot the CPA DOS, 
the ATA DOS, and the numerically calculated DOS, with 
the left Y-axis representing their values. We add the CPA 
applicability criterion by plotting the small parameter of 
the series R(e) (solid line), with the right Y-axis show- 
ing its values in the same figures. We also designated by 
the triangle the Fermi level position, obtained from the 
numerical data for the impure system. 

For the low concentrations (i.e. those that are less 
than c r ), analytical curves, namely the CPA DOS and 
the ATA DOS, perfectly fit the numerical histogram. The 
DOS only slightly deviates from the conventional Dirac 
DOS mainly because of the shift towards positive en- 



ergies. The applicability criterion R(e) < 0.5 is satisfied 
practically at all energies within the chosen window, R(e) 
is small and characteristically contains two maxima. The 
sharp one corresponds to the van Hove singularity and 
predicts failure of analytical approximations in the Dirac 
point vicinity. Less sharp one is due to the a 2 factor. 

When the impurity concentration is increased approx- 
imately to the critical value c r , maxima of small param- 
eter R(e) show the tendency to merge together into a 
single maximum, which height goes beyond the 0.5 value 
(as it was shown in Ref. 9). This event indicates the on- 
set of the spectrum rearrangement, providing a reference 
point for the critical concentration at the given Vl- In 
the domain with heightened values of R(e) the discrep- 
ancy between the CPA DOS and the simulation results 
is more clearly expressed. 

Figures also show that the perturbation Vl — 4t is 
marginal as the resonance state appears at the periphery 
of the region, where the Dirac approximation ^ works 



well. Even for low concentrations some divergence is seen 
between analytical and numerical curves at the edges 
of energy domain considered. Impurity resonances are 
smeared out and, therefore, cannot be readily discerned 
for a perturbation of this strength. Likewise, for such a 
Vl the impurity resonance is not well defined in the single 
impurity LDOS. It should be mentioned that in the weak 
scattering regime (Vl < W) the spectrum rearrangement 
process does not take place at all (as it was evident in 
Ref.0, when the average DOS maintained the mere rigid 
shift with increasing impurity concentration). 

With an increase in Vl the impurity resonance be- 
comes well defined. It has been obvious beforehand that 
the CPA DOS should not contain any sharp features in 
a contrast to the ATA DOS. Because of the absence of 
self-consistency the ATA gravitates more to the single- 
impurity resonance. It is clearly seen in Fig. [4] and Fig. [5] 
that the ATA DOS quite correctly reproduces the reso- 
nance peak at the impurity concentrations that are close 
to the critical one. The larger is the impurity potential 
Vl, the better the ATA curve fits numerical data. 

This coincidence, however, is not pertained to the im- 
purity concentrations that exceed the critical concentra- 
tion of the spectrum rearrangement. With an increase in 
the impurity concentration the maximum position shifts 
in the direction of negative energies from the energy of 
the single-impurity resonance and the second, accom- 
panying maximum is coming forth. While the shift of 
the primary maximum from the Lifshitz equation root 
is considerable, the maximum in the ATA DOS remains 
still at the single- impurity resonance energy. The afore- 
mentioned irregular structure in the DOS is the most 
intriguing feature of the spectrum which cannot now be 
interpreted with the help of the available CPA or ATA 
approximations. The maxima in the DOS are located 
within a large domain, in which the CPA DOS does not 
follow the simulated curve and, correspondingly, the CPA 
validity criterion is not satisfied (R(e) > 0.5). This do- 
main covers the single-impurity resonance and the mini- 
mum in the DOS, at which valence and conduction bands 
are docking each other. 

In addition, to study the character of states within this 
domain, we calculated the inverse participation ratio 7 ' 
P(e) as a localization criterion in a system of a smaller 
size: 



x 10 



(11) 



where summation runs over all lattice sites. Even though 
the comparison for systems of different size is not in- 
cluded in the current article, we should mention that the 
hight of P{e) curves does not diminish with increasing 
the size of the system suggesting electron localization. 
The dynamics of the inverse participation ratio with in- 
creasing the impurity concentration is given in Fig. [6] 
Here again periodical boundary conditions were used at 
zigzag edges of the sample to get rid of the sharp peak at 
e = 0, which is related to boundary states. Chosen con- 





0.015 

CL 

0.005 




0.015 



c = 0.00075 





c = 0.003 


.1 111 IL.ii 


dilU..,,, . 




-0.25 -0.2 



-0.1 -0.05 



FIG. 6: (color online) Inverse participation ratio for Vl 
and different impurity concentrations. 
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FIG. 7: (color online) Fragments of eigenstates for Vl = 
and c = 0.012 are depicted. These configurations of impuriti 
are characteristic for the second peak in the DOS. 
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centrations repeat those from Fig. [4] that corresponds to 
the same impurity potential. A radical localization inten- 
sification after the spectrum rearrangement is evidently 
seen. States arc showing a tendency for their localization 
in the very region, in which the CPA is not valid. This 
fact also confirms a close connection between the CPA ap- 
plicability criterion and the Ioffe-Regel criterion^ When 
comparing Fig. [5] to Fig. [|J in particular at c = 0.012, it 
is obvious that the largest values of P(e) match the peaks 
in the DOS. 

A pair of states corresponding to the sharp peaks in 
the inverse participation ratio graph are shown in Figs. [7] 
and [8] The states corresponding to the first (counting 
from e = 0) maximum in the DOS at e w — 0.14i are 
mostly represented by relatively distant impurity pairs 
and triads. Equally challenging is the origin of the second 
DOS peak at e w — 0.19t. The visualization of the wave- 
function belonging to this region is provided in Fig. [7] 
It shows that this peak is largely due to the characteris- 
tic pattern of impurities, which is depicted in the same 
figure. It is worth mentioning that impurity atoms are lo- 
cated on one sublattice for these strongly localized states, 
while the ^-function is concentrated on the other sublat- 
tice. It resembles the situation with a double impurity 3 
and can be attributed to the relation |.go Q ™J "C |go Q Ti^| 
for \e\ — ► 0. To summarize, when the critical concen- 
tration of impurities is exceeded, a quasigap filled with 
localized states is developing in the graphene's spectrum 
because of the impurity cluster effects. Since the CPA 
does not account for cluster effects and scatterings by im- 
purity clusters dominate within this quasigap, the CPA 
is not applicable in this region. 

Numerical results show that Dirac point as such is 
eliminated from the spectrum when the critical concen- 
tration c r is reached. Consequently, it is not justifiable to 
speak about the Dirac point existence, albeit the impu- 
rity concentration can be relatively low (as low as c r is) . 
The CPA and the ATA do not correctly describe the DOS 
minimum between the valence band and the conduction 
band for c > c r . Normally, the Dirac point coincides with 
the Fermi energy for the pure (or undoped) graphenc. 
However, at the finite concentration of impurities situa- 
tion changes drastically. Dean's approach allows to track 
the position of the Fermi level, since it outputs the total 
number of states located below any given energy. The 
Fermi level monitoring revealed that its position shifts in 
the positive direction away from the DOS minimum. The 
greater is the impurity concentration the more prominent 
is the Fermi energy shift from the DOS minimum, which 
transforms graphenc into a "doped" conductor without 
any gate voltage applied. Should impurity atoms bring 
additional electrons in the system, the doping effect will 
be increased. 

This shift can be explained by the following uncompli- 
cated consideration. Without the conduction band, the 
valence band should develop a strong tail above its up- 
per edge for a positive impurity potential Vl > 0. Since 
the conduction band is not separated from the valence 
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FIG. 8: (color online) Fragments of the eigenstates for Vl = St 
and c = 0.012 are depicted. They correspond to the first 
(closest to the Dirac point) peak in the DOS. 



band, this tail falls inside the conduction band produc- 
ing excessive states above the DOS minimum. Keep- 
ing in mind that total number of states within the sin- 
gle band should be preserved because of the sum rules, 
this expansion of the valence band into the conduction 
band yields a deficit of states below the DOS minimum. 
At the constant amount of carriers, extra electrons are 
flowing out to the conduction band altering the Fermi 
level position. With increasing impurity concentration 
the tail is gradually becoming more pronounced, which 
consequently makes the shift in the Fermi energy more 
apparent. 



IV. CONCLUSION 

A comprehensive analysis of the spectrum rearrange- 
ment in graphene with substitutional impurities by nu- 
merical simulation has been carried out. We studied the 
DOS near the Fermi level in graphene for a set of impu- 
rity potentials and impurity concentrations. 
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It was demonstrated that indeed a certain characteris- 
tic concentration of impurities can be specified, at which 
the graphene's spectrum undergoes a qualitative change. 
This critical impurity concentration is associated with 
the spatial overlap of individual impurity states. In a 
turn, it has been established that the cardinal modifi- 
cation of the spectrum is manifested by the opening of 
the filled with highly localized states quasigap around 
the impurity resonance energy. The cluster effects were 
found to be responsible for the quasigap formation. Pair 
impurity states representing the most prominent peaks 
in the DOS within the quasigap have been visualized, 
which emphasized the dominance of scatterings on im- 
purity clusters inside the quasigap. Aforesaid confirmed 
the predicted scenario of the spectrum rearrangement in 
graphene. 

A comparison of the CPA DOS with the numerically 
simulated DOS supported the suggested CPA applica- 
bility criterion and its efficacy as an instrument in the 
description of the spectrum rearrangement passage. As 
well, intimate correlation between the CPA validity and 



the degree of electron localization has been revealed. 
That is, inside the quasigap, in which cluster effects are 
essential and states are localized, the CPA is not reliable. 

Furthermore, we report about the phenomenon of the 
Fermi level shift from the DOS minimum — a kind of 
a self-doping, which alters the conductivity of impure 
graphene without gate voltage variation even in the case 
of neutral impurities. 
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